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Abstract. A stereo correlation method on the object domain is proposed to generate 
the accurate and dense Digital Elevation Models (DEMs) from lunar orbital imagery. 
The NASA Ames Intelligent Robotics Group (IRG) aims to produce high-quality 
terrain reconstructions of the Moon from Apollo Metric Camera (AMC) data. In 
particular, IRG makes use of a stereo vision process, the Ames Stereo Pipeline (ASP), 
to automatically generate DEMs from consecutive AMC image pairs. Given camera 
parameters of an image pair from bundle adjustment in ASP, a correlation window is 
defined on the terrain with the predefined surface normal of a post rather than image 
domain. The squared error of back-projected images on the local terrain is minimized 
with respect to the post elevation. This single dimensional optimization is solved 
efficiently and improves the accuracy of the elevation estimate. 


1. Introduction 

Topographical maps are an essential tool for scientists interested in exploring and learning 
more about planetary bodies like the moon or mars. These maps allow scientists do 
everything from identifying geological phenomena to identifying potential landing sites for 
probes or spacecraft. 

Satellites and other spacecraft that visit planetary bodies of interest are usually equipped 
with a variety of sensors, some of which can be used to recover the topography of the 
planetary surface. LiDAR (Light Detection And Ranging) sensors give sparse (but highly 
accurate) measurements at periodic points called “posts.” Raw images that are captured as 
the satellite orbits the planetary body can also be processed to create highly detailed 
topographical maps. By registering and aligning these two data sets, maps can be created 
that are both dense and accurate. 

Given two images of the same scene taken from slightly different perspectives, the 
relative shift of objects from frame to frame (known as “disparity”) is related to distance of 
the object: far objects appear to move less than close objects. This phenomenon should be 
familiar, since the human brain uses this relationship to create depth perception from the 
differences in perspective seen by both eyes. Similarly, depth information from orbital 
imagery can be recovered in areas where the images overlap by matching points in the 
images that correspond to the same 3D location and measuring their disparity. By using the 




(a) The Mapping Cameras System (b) Apollo 15 Orbit 33 

Figure 1: AMC Data System, (a) The AMC captures a series of pictures of the Moon's surface, (b) 
Satellite station positions for Apollo Orbit 33 visualized in Google Moon. 


disparity along with the location and orientation of the spacecraft as well as a mathematical 
model for the lens of the camera, the precise location of the 3D point can be calculated. This 
technique of using matching points between images to calculate 3D locations is known as 
“stereophotogrammetry” [6]. 

Before computers automated this task, points in orbital images were manually matched 
using mechanical stereoplotters. Now, computers can perform this once long and arduous 
process with minimal human interaction. For example, the Ames Stereo Pipeline [2, 5] is a 
collection of cartographic and stereogrammetric tools for automatically producing digital 
elevation maps (DEMs) from orbital images acquired with the Apollo Metric Camera 
(AMC) during Apollo 15-17 (Figure 1). 

A new stereo correlation method is proposed on the terrain model to improve the 
accuracy of the DEM. This paper will address this problem by finding robust elevation from 
an image pair that minimizes the discrepancy between two back-projected patches on the 
local planar terrain. The proposed method determines the elevation by minimizing the 
squared error of the back-projected image patches. The accuracy of DEMs produced by 
IRG will thus be improved. The reconstructed DEMs from lunar orbital imagery are 
presented. 


2. Ames Stereo Pipeline 



Figure 2: Dataflow of the Ames Stereo 
Pipeline. Preprocessing includes the 
registration and filtering of the image pair. 
A stereo correlator (disparity initialization 
and sub-pixel refinement) constructs the 
disparity map based on normalized cross 
correlation. DEMs are generated by a 
triangulation method in which corrected 
camera poses are used by bundle 
adjustment. 


The Ames Stereo Pipeline (ASP) is the stereogrammetric platform that was designed to 
process stereo imagery captured by NASA spacecraft and produce cartographic products 
since the majority of the AMC images have stereo companions [10] [18]. The entire stereo 
correlation process, from an image pair to DEM, can be viewed as a multistage pipeline 
(Figure 2). At the first step, preprocessing includes the registration to align image pairs and 
filtering to enhance the images for better matching. Triangulation is used at the last step to 
generate a DEM from the correspondences. 


2.1. Bundle Adjustment 

Before stereo correlation is performed, Bundle Adjustment (BA) corrects the 
three-dimensional postures of cameras and the locations of the objects simultaneously to 
minimize the error between the estimated location of the objects and their actual location in 
the images. Camera position and orientation errors affect the accuracy of DEMs produced 
by the ASP. If they are not corrected, these uncertainties will result in systematic errors in 
the overall position and slope of the DEMs. BA ensures that observations in multiple 
different images of a single ground feature are self-consistent. 

In BA the position and orientation of each camera station are determined jointly with the 
3D position of a set of image tie-points points chosen in the overlapping regions between 
images. Tie-points are automatically extracted using the SURE robust feature extraction 
algorithm [14]. Outliers are rejected using the random sample consensus (RANSAC) 
method [15]. The BA in ASP determines the best camera parameters that minimize the 
reprojection error [16]. The Levenberg-Marquardt algorithm is used to optimize the cost 
funciton [17]. 



(d) Parabola Refined Map (e) Cauchy Refined Map (f) Bayesian Refined Map 

Figure 3: Stereo Correlation, (a-b) An image pair, (c-f) Horizontal disparity maps, (c) The fast 
discrete correlator constructs the coarse disparity map. (d-f) Refined disparity maps from the 
initialized map. (f) Bayesian sub-pixel correlator generates a smoother map than the others. 


2.2. Stereo Correlation 

Stereo correlation, which is the process at the heart of ASP, computes pixel 
correspondences of the image pair (Figure 3). The map of these correspondences is called a 
disparity map. The best match is determined by applying a cost function that compares the 
two windows in the image pair. The normalized cross correlation is robust to slight lighting 
and contrast variation in between a pair of images [11]. For large images, this is 
computationally very expensive, so the correlation process is split into two stages. (1) The 
disparity map initialization step computes coarse correspondences using a multi- scale 





(a) Left Image (b) Right Image (c) Rendered DEM 

Figure 4: Generation of DEM. (a) and (b) Apollo Metric Camera image pair of Apollo 15 site, (c) A 
DEM of Hadley Rille is rendered from an image pair. 


search that is highly optimized for speed (Figure 3c). (2) Correlation itself is carried out by 
sliding a small, square template window from the left image over the specified search 
region of the right image (Figure 3d-f). 

Refining the initialized disparity map to sub-pixel accuracy is crucial and necessary for 
processing real-world data sets [13]. The Bayesian expectation maximization (EM) 
weighted affine adaptive window correlator was developed to produce high quality stereo 
matches that exhibit a high degree of immunity to image noise (Figure 3f). The Bayesian 
EM sub-pixel correlator also features a deformable template window that can be rotated, 
scaled, and translated as it zeros in on the correct match. This affine-adaptive window is 
essential for computing accurate matches on crater or canyon walls, and on other areas with 
significant perspective distortion due to foreshortening. A Bayesian model that treats the 
parameters as random variables was developed in an EM framework [21]. This statistical 
model includes a Gaussian mixture component to model image noise that is the basis for the 
robustness of the algorithm. The resulting DEM is obtained by triangulation method 
(Figure 4). 

3. Orthographic Stereo Correlation 


To improve the DEMs produced by IRG, we proposed to define error measure in object 
domain rather than image domain for lunar orbital images. We propose to use a planar 
approximation of lunar terrain, which define the relationship of two views by a homography 
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(c) Lunar reflectance (d) Linear reflectance 

Figure 5: Linear Approximation of Lunar Terrain and Reflectance, (a) Smooth local 
terrain in a small field of view is approximated by (b) planar one. (c) The lunar 
reflectance model is approximated by a linear function, (d) The linear reflectance is a 
linear function of surface albedo (A). 


representation. To model lunar reflectance, we propose to use a linear reflectance 
approximation. The statistical behavior of the photons is model by the normal distribution 
to derive the cost function that compares the two view windows. The proposed method will 
replace the pair-wise sub-pixel refinement and triangulation currently used in ASP. 

A linear approximation of the lunar terrain and reflectance simplifies the correlation 
function. The lunar surface is smooth, but not flat (Figure 5a) and has its own reflectance 
(Figure 5c). The terrain is assumed to be locally planar because the correlation patches are 
determined to be small enough to cover the planar region (Figure 5b). Similarly, the lunar 
reflectance is assumed to be locally linear because geometric and photometric changes are 
small enough to form a linear relationship (Figure 5d). 

A planar approximation of the terrain provides a homography representation of two-view 
correspondences. Suppose we have a textured plane that approximates the local terrain 
(Figure 7). Let n be the surface plane with the normal vector n and distance d from the 
origin, i.e., n \ vl x- d . Let I be the orbital image and / be the orthographic image of 
/ onto n . The surface albedo having linear relationship with / , i.e., /( z) = b + cl( z) 


Figure 7: Homography Representation of Stere 
Correspondences. Homographies of the image 
planes induced by the local plane /r, 
H,:/; I — ^ 1 . for i — 1,2 , determine 

correspondences of the image pair. 


for all t E: 7i . Suppose we have an image pair viewing the same scene on n with projection 
matrices P 1 and P 2 . 

To determine the geometric and reflectance parameters, we propose to minimize the 
squared error of all corresponding patches on the terrain. From the fact that the quantity of 
incident photons from the scene follows the normal distribution, the negative log-likelihood 
of scene radiance is derived: 

L{ b, c, d; t) = J g(z — t) {/ (z) - f 2 (z)} 2 dz . (1) 



(a) AS15-M-1090 


(b) AS15-M-1091 


Figure 6: A Stereo Image Pair. Two images from Apollo 15 Orbit 33 imagery have different 
illumination and noises such as dust and missing data. 



where f.(z) = b i +c i I.(z) , b and c are coefficient vectors of linear reflectance, and 
(•) T =g( t)** is the local Gaussian convolution operator. The optimizer of (1) may not 
exist in a closed form, but some parameters are simply optimized from the other parameters. 
To minimize the squared error, we propose to employ an alternating optimization scheme. 
First, we solve for b and c because there is a closed- form solution given d in (1). Then 
we solve for d in (1), while keeping surface albedo constant with the solutions from 
(l).The initialized disparity maps from ASP can be used to calculate the initial value of d . 

4. Experimental Results 

The orthographic stereo correlator is implemented based on the NASA Vision Workbench 
(VW). The NASA VW is a general purpose image processing and computer vision library 
developed by the IRG at the NASA Ames Research Center. The linear reflectance 
coefficients are estimated in the nested optimization process. Since this is a single 
dimensional optimization with respect to the elevation, we adopted a hybrid golden section 
and parabolic interpolation method. Figure 6 shows an image pair from Apollo metric 
images. A Gaussian window with scaling factor 5 is used for correlation window by 33*33 
pixels. Back-projected images are generated from an image pair with bicubic interpolation. 
The surface normal is determined by the radial direction of the local terrain because the 
lunar terrain is smooth enough. 

Figure 8a shows a stereo DEM constructed from integer disparity map. We can observe the 
pixel locking effect in the initial DEM with many holes. With this initial DEM, a refined 
DEM is reconstructed in Figure 8b. As you can see in the figure, the refined DEM gets rid 
of pixel locking artifact completely. There are small bumps and holes because the initial 
DEM has the large error and optimizer falls in local minima. Even though the surface 
normal is assumed roughly to have the radial direction of each post, the reconstructed 
terrain is reconstructed smooth enough and robust. 

5. Conclusion 

An orthographic stereo correlation method was proposed to reconstruct the accurate and 
dense Digital Elevation Models (DEMs) from lunar orbital imagery. The proposed method 
addresses this problem by making use of linearity in geometry and photometry to improve 
both the accuracy and robustness of the stereo correlation process. Given camera 
parameters of an image pair and an initial DEM, the DEM is refined to acquire the sub-pixel 
accuracy. A correlation window on the terrain with the predefined surface normal of a post 
is used to define the squared error of back-projected images on the local terrain. The 
squared error is minimized with respect to the post elevation. This forms a single 
dimensional optimization by nesting the linear reflectance optimization. For the future 
work, the surface normal can be estimated together with the elevation. 




Figure 8: Dense DEM Reconstruction. Color-mapped, hill- shaded DEM was reconstructed from 
from an image pair in Figure 6. 


6. References 

1. Noble, S.K., et al., The Lunar Mapping and Modeling Project. LPI Contributions, 2009. 1515: pp. 
48. 


2. Mordohai, P. and G. Medioni. Dense multiple view stereo with general camera placement using 
tensor voting, in IEEE 3DPVT. 2004: pp. 725-732. 

3. Thompson, C.M., Robust photo-topography by fusing shape-from-s hading and stereo. 1992, 
Massachusetts Institute of Technology. 

4. Price, K., Annotated computer vision bibliography. 1995: Institute for Robotic and Intelligent 
System School of Engineering (IRIS), University of Southern California. 

5. Heipke, C. and C. Piechullek. Toward surface reconstruction using multi-image shape from 
shading, in ISPRS. 1994: pp. 361-369. 

6. Cryer, J.E., P.S. Tsai, and M. Shah, Integration of shape from shading and stereo. Pattern 
recognition, 1995. 28(7): pp. 1033-1043. 

7. Hapke, B., Bidirectional reflectance spectroscopy, 1, Theory. J. Geophys. Res, 1981. 86 : pp. 
3039-3054. 

8. Minnaert, M., Photometry of the Moon. Planets and Satellites, 1961. pp. 213. 

9. Buratti, B. and J. Veverka, Voyager photometry of Rhea, Dione, Tethys, Enceladus and Mimas. 
Icarus, 1984. 58(2): pp. 254-264. 

10. Broxton, M.J., et al., The Ames Stereo Pipeline: NASA’s Open Source Automated 
Stereo grammetry Software. 2009, NASA Ames Research Center. 

11. Menard, C., Robust Stereo and Adaptive Matching in Correlation Scale-Space. 1997, Institute 
of Automation, Vienna Institute of Technology. 

12. Sun, C., Fast stereo matching using rectangular subregioning and 3D maximum-surface 
techniques. International Journal of Computer Vision, 2002. 47(1): pp. 99-117. 

1 3 . Nefian, A., et al. A Bayesian Formulation for Subpixel Refinement in Stereo Orbital Imagery, in 
International Conference on Image Processing. 2009. Cairo, Egypt. 

14. Bay, H., et al., Speeded-up robust features (SURF). Computer Vision and Image Understanding, 
2008. 110 ( 3 ): pp. 346-359. 

15. Fischler, M.A. and R.C. Bolles, Random sample consensus: A paradigm for model fitting with 
applications to image analysis and automated cartography. Commun. Assoc. Comp. Mach,, 1981. 
24(6): pp. 381-395. 

16. Triggs, B., et al., Bundle adjustment - a modern synthesis. Lecture Notes in Computer Science, 
1999. pp. 298-372. 

17. Hartley, R. and A. Zisserman, Multiple view geometry in computer vision. 2003: Cambridge 
University Press. 

18. Taemin Kim, Zachary Moratto and Ara V. Nefian, Robust Mosaicking of Stereo Digital 
Elevation Models from the Ames Stereo Pipeline, 

19. F. Pedersini, P. Pigazzini, A. Sard, S. Tubaro: 3D Area Matching with Arbitrary Multiview 
Geometry. EURASIP Signal Processing: Image Communication - Special Issue on 3D Video 
Technology", Elsevier, Vol. 14, No 1-2, October 1998, pp.71-94 

20. M. J. Broxton, A. V. Nefian, Z. Moratto, T. Kim, M. Lundy, and A. V. Segal. 3D Lunar Terrain 
Reconstruction from Apollo Images. In Proceedings of the 5th IS VC 2009 

21. A. Nefian, K. Husmann, M. Broxton, V. To, M. Lundy, M. Hancher. A Bayesian Formulation 
for Sub-pixel Refinement in Stereo Orbital Imagery, in Proceedings of the 2009 IEEE ICIP 2009 



